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We explore observational constraints on possible deviations from Newtonian gravity by means of 
large-scale clustering of galaxies. We measure the power spectrum and the bispectrum of Sloan Dig- 
ital Sky Survey galaxies and compare the result with predictions in an empirical model of modified 
gravity. Our model assumes an additional Yukawa-like term with two parameters that characterize 
the amplitude and the length scale of the modified gravity. The model predictions are calculated 
using two methods; the second-order perturbation theory and direct A'^-body simulations. These 
methods allow us to study non- linear evolution of large-scale structure. Using the simulation re- 
sults, we find that perturbation theory provides reliable estimates for the power spectrum and the 
bispectrum in the modified Newtonian model. We also construct mock galaxy catalogues from 
the simulations, and derive constraints on the amplitude and the length scale of deviations from 
Newtonian gravity. The resulting constraints from power spectrum are consistent with those ob- 
tained in our earlier work, indicating the validity of the previous empirical modeling of gravitational 
nonlinearity in the modified Newtonian model. If linear biasing is adopted, the bispectrum of the 
SDSS galaxies yields constraints very similar to those from the power spectrum. If we allow for the 
nonlinear biasing instead, we find that the ratio of the quadratic to linear biasing coefficients, 62/61, 
should satisfy —0.4 < 62/61 < 0.3 in the modified Newtonian model. 

PACS numbers: 04.50.-|-h 98.65.-r 98.80.Es 



I. INTRODUCTION 

Recent measurements of cosmic microwave background 
radiation angular power spectrum [l|, d, 0] strongly sup- 
port the "standard" cosmological model in which the en- 
ergy content of the universe is dominated by dark energy 
(very close to Einstein's cosmological constant, A) and 
cold dark matter (CDM). Such ACDM universes are also 
in good agreement with independent datasets ofgalaxy 
clustering d, and distant Type la supernovae [a 0, 111 ■ 
Thus the basic framework for the theory of structure for- 
mation in the universe is firmly established. Nevertheless 
the nature and the physical origin of dark energy remain 
to be understood. 

The apparent accelerating expansion of the universe is 
conventionally interpreted in terms of a source of repul- 
sive force (dark energy), but can be explained by mod- 
ifying Newton's law of gravity on cosmological scales as 
well. The latter resolution has been seriously considered 
recently. For example, Dvali, Gabadadze and Porrati 
(DGP) 0, ^ propose that gravity leaking into extra 
dimensions drives the observed accelerating expansion. 
Other such models include modified Newtonian dynam- 
ics (MOND) [ll|,[ll,[ll and ghost condensation [3,111. 
Intriguingly, all of these alternative models predict some 
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deviation from conventional Newtonian gravity at cos- 
mological scales. 

Indeed, while the validity of Newtonian gravity is 
tested to high precision up to the scale of the solar sys- 
tem (~ lO^'^ m), there have been no rigorous tests at 
sub-millimeter and over scales beyond the solar system 
[3, 113, Ell • It has been suggested that large-scale galaxy 
clustering can be used to constrain non-Newtonian mod- 
els of gravity in principle, but it became feasible only 
recently with accurate measurements of galaxy clustering 
in large rcdshift surveys [1, Q . 

In our earlier work po| (Paper I), we put quantitative 
constraints on deviations from Newtonian gravity at cos- 
mological scales under the assumption that the deviation 
can be described in a simple parametric form; we adopted 
an empirical Yukawa-like term for the modified gravity, 
and calculated the galaxy-galaxy power spectrum semi- 
analytically. (See also Ref. [21| for similar argument.) 
By comparing thepredicted power spectrum with that 
of SDSS galaxies [1], we derived quantitative, although 
still conditional, constraints on deviations from Newton's 
law of gravity. 

In this paper, we improve our previous work by 
performing non-linear cosmological simulations and by 
exploiting a higher-order statistic, bispectrum. Since 
bispectrum is sensitive to clustering in the non-linear 
regime, it is expected to provide complementary con- 
straints at mega-parsec scales to that obtained from 
power spectrum analysis. We use direct iV-body sim- 
ulations to test the accuracy of our semi-analytic calcu- 
lations and to reinforce our conclusions. 
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The rest of the paper is organized as follows. Our 
model assumptions are described in Sec. [Ill We derive 
power spectrum and bispectrum from perturbation the- 
ory in modified Newtonian model in Sec. IIIII and HVl We 
perform A^-body simulations and construct mock samples 
of volume-limited SDSS galaxies for direct comparison 
with the observational data. Details of the simulations 
are described in Sec. |Vl The results of perturbation the- 
ory and the simulations are discussed in Sec. IVII and lVIII 
Finally Sec. IVIIII concludes the present analysis. 



where Sk g^iaxy ^^'^ ^fc fractional fluctuation of galaxy 
number and mass density, bi and 62 are linear and 
quadratic biasing parameters. We consider only linear 
bias (i.e., &2 = 0) when we use power spectrum, whereas 
we consider both bi and 62 for analyses using bispectrum. 
To derive constraints on a and A, bi is treated as a free 
parameter to adjust the overall clustering amplitude. 

III. POWER SPECTRUM ANALYSIS 



II. MODEL ASSUMPTIONS 

In this section, we briefly summarize our model and 
a set of assumptions. Further details may be found in 
Paper I. 

We consider a modified Newtonian model for which 
gravitational potential is given by 
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where Gn denotes (conventional) Newton's constant of 
gravity. The above model corresponds to Model II in Pa- 
per I, on which we focus throughout the following anal- 
ysis. The deviation from the Newtonian gravity in this 
model is characterized by two parameters, a and A; a is 
the dimensionless amplitude of the deviation and A is the 
characteristic length scale. Note that A is defined in the 
proper length, rather than in the comoving length. 

It is important to note that, although we consider de- 
viations from Newtonian gravity at mega-parsec scales, 
we still assume that the global cosmic expansion is unaf- 
fected by the deviations. Namely, we assume that gen- 
eral relativity is valid on horizon scales and thus the cos- 
mic expansion is described by the standard Friedmann 
equation. Strictly speaking, these two assumptions may 
be in conflict with modified gravity models in general 
[H, m, [13, m, nil . To account for the existing data such 
as SNela and CMB, however, the cosmic expansion law 
can hardly change in practice. This is why we adopt the 
conventional Friedmann equation even in this analysis. 
For the same reason, we use conventional matter trans- 
fer function as initial condition of dark matter adopting 
the background cosmology defined by the standard set of 
cosmological parameters, f2ni=0.3, r2b=0.04, f^A = 0.7, 
and the Hubble constant at present h = 0.7 in units of 
100km s^^ Mpc^^. See Paper I for further discussion on 
this point. 

In order to make a direct comparison between the clus- 
tering of SDSS galaxies and our model predictions, we 
need to assume a biasing relation for the distribution of 
galaxies and that of matter. For this purpose, we adopt 
a commonly adopted deterministic relation: 
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In Fourier space, the modified gravitational potential 
in Eq. ((T|) can be written as 



[A,$(a;)]^ = 47rGNa'p 
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where x is in the comoving coordinate, k is the comoving 
wave-number, and a is the scale factor normalized unity 
at the present epoch. 

For the potential of Eq. ([S]), the evolution equation for 
density perturbations is written as 



with 
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where H{a) is the Hubble parameter, and 5^^^^ denotes the 
linear term in density fluctuations [see Eq. (|16p below]. 
Note that even the linear perturbation equation becomes 
dependent on k in the modified gravity model. 

Next, the linear power spectrum Pi_,{k) at present is 
given by 



PL(fc;a,A) = Ar^(fc)fc" 6'^^\a = l]a,\) 



(7) 



where T{k) is the matter transfer function, and n is the 
spectral index of the primordial power spectrum which 
we set to be unity. We use the fitting formula of Eisen- 
stcin and Hu [23| for T{k). It should be emphasized 
here that we fix the amplitude A so that the rms value of 
the top-hat mass fluctuations at 8ft.~^Mpc, cts, equals 0.9 
when a = and n = 1. The actual value of a% in our 
modified gravity model may be slightly different because 

of the factor 5l,'{a = l;a. A) in Eq. ([7]). However, 

the difference in the overall amplitude is unimportant 
because we have an additional freedom to adjust the pre- 
dicted amplitude via the biasing relation [Eq. ([2|)]. 
In Paper I we used the Peacock-Dodds prescription 
to convert the linear power spectrum to nonlinear 
one. It turned out that in doing so we used incorrectly 
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the growth factor g{^) and the tilt of Unear power spec- 
trum ni^{ki^) given in the case of Newtonian models. We 
made sure later that the above mistake did not change 
the final power spectra very much as long as the Peacock- 
Dodds prescription is valid. In the present paper, we 
also confirm the validity of the Peacock-Dodds approach 
in non-Newtonian models using A^-body simulations di- 
rectly (see Sec. EH). 



These equations can be solved recursively. Let us first 
decompose Sk and perturbatively, 



5f +4^^ 



Ok = Ofc 



(16) 
(17) 



Differentiating Eq. (|11[1 and substituting to Eq. p4[l to 
eliminate Uk, we obtain 



IV. PERTURBATION THEORY AND 
BISPECTRUM 

In this section, we describe the second order perturba- 
tion theory and its application to bispectrum. The earlier 
formulation of cosmological perturbation in the Newto- 
nian model may be found in [2^, [s^, [3l|, [11] . Bernardeau 
[33I I developed a formulation of second order perturbation 
theory in non-Newtonian models. We apply the method 
to the modified potential in Eq. 

The basic equations are given by 
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where the over-dot denotes the derivative with respect 
to time, «'(= ai*) is the peculiar velocity, and <& is the 
gravitational potential. We define velocity divergence: 
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where A*-"^ and i?*-"^ are the source terms of the nth- 
order: 
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Consider first the lowest order, n = 1. Since Ad) = 
= 0, Eq. IIHI) reduces to 
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0, 
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Equations 
aSk + Ufc = - 



and 



uix,t)=dy{x,t). (10) 
and ^ in Fourier space reduce to 
1 
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which is equivalent to Eq. We denote the growing 



mode of the solution of Eq. (gT]) by Dl'^t). Note that, 
in non- Newtonian models, the solution D^\t) is gener- 
ally dependent on scale k, in contrast to the conventional 
Newtonian case. The linear solution 5^^' is given by 



5i'^^D['^5Uk), 



(22) 



where (5ini(fc) is the initial fractional density. 

The corresponding linear solution for Uk is obtained 
from Eq. PT|) as 



(23) 



Solutions at the next order, n — 2, are more compli- 
cated. Eq. for 71 = 2 is written explicitly as 
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= / d^pd^qS''ip + q-k){6MSini{q)[So{p,q,t)Poi^i) + Slip,q,t)Vl{^l^ (24) 

Soip,q,t) = + ^) 4^^^^^^ + IdI^'D^,''^ (25) 

^ib^'^^^) = + + + ^) ^i'^^'.'^^ (26) 

52(p,g,t) = ^I)(i)l)(i), (27) 



where 6^{k) is the Delta function and P/(/x) are the Leg- 
endre polynomials: 



Po(/^) = l, Pi(A*)=/i, P2(aO = ^ (3m' - 1) , (28) 



Equation 



pq 
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pq 

has an implicit solution of the form: 
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y" d^p£q6^{p + q-k) [(5i„i(p)(5ini(9) 



{ro(p, g, t)VQ{^l) + Tlip, q, t)Vi{ti) + T^ip, q, t)V2if^)}] , 

(30) 

where the functions Ti(p,q,t) satisfy 

V\p+g\T,ip,q,t)^S,ip,q,t) fori = 0,1, 2. (31) 

We note that expressions for the second-order solutions 
given in [3^ contain some typographical errors which are 
corrected in our above expressions. 

These results enable us to compute the bispectrum in 
the leading order. The bispectrum is defined as 

(^(fel)^(fc2)^(fc3)) = (27r)^B(fci,fc2,fc3)<5^(fci+fc2 + fc3). 

(32) 

The leading-order terms of the left-hand-side of the above 
equation are given by 

{5iki)Sik2)Sik3)) = {S^^\k,)S^^\k2)S^'\k3)) 

+ cj/c.(l,2,3). (33) 
Therefore the bispectrum reduces to 
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where Phu{k) = (|(5ini(fc)|^). In what follows, we write the 
bispectrum simply as B{ki, k2) adopting the condition of 
fcg = -fc2 -fci [Eq. (pH)]. 

To compute the bispectrum, we solve Eq. (j3ip nu- 
merically for each pair of (p, q), together with the linear 
perturbation equation (j2ip . At sufficiently early epochs 



{zi ^ 1), D^j^\zi) is simply given by the growth rate in 
the Newtonian case (see Paper I) . Similarly, Ti are given 
by 



To(p,q,z,) = 11(1 + z,)-^ 
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Ti(p,9,z,) = l + (l + z,)-2, (36) 

2 \p qj 
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SIMULATION AND OBSERVATIONAL 
DATA 



A. N-body Simulations 

We use the cosmological A^-body solver TPM-1.1 
[131 in its PM-only mode. We run six realizations 
each for simulation box-sizes of ibox = 500/i~^Mpc, 
and 1000/i~^Mpc with the following parameters: a = 
-1.0,-0.8,-0.5,-0.2,0.0,0.2,0.5,0.8, and 1.0, A = 
2, 5, 8, 10, 12, 15, 20, and ZQh-^M^c. We use the fitting 
formula for the matter transfer function, equation (28) 
^ (31) of the ref. [l^l, that ignores the baryon acoustic 
oscillation effect. We start the simulations at z = 50. All 
the simulations employ Np = 128^ particles. 

To simulate structure formation in the non-Newtonian 
model, we need to modify the Green function of the 
Laplacian, Q. For a density field p defined on a three- 
dimensional wave- number grid (p, g, r), the gravitational 
potential in real space is evaluated to be 

(^(/, m, n) — ^ Qp.q.rPp.q.r cxp [2ni{pl + qm + rn)/M] , 

p^q^r—O 

where /, m, n are position integers in real space with M 
being the number of grids per dimension (we follow the 
notation in Efstathiou et al. (ssj). 

The Green function in the original TPM code that as- 
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sumcs the conventional Newtonian gravity is given by 

{0, I = in = n = 0; 

-7T/ {M2 [sin^iTTp/M) 
+ sin^(7rg/M) + sin^(7rr/A/)] } otherwise; 

(39) 

which is derived from the seven-point finite-difference ap- 
proximation. 

Taking account of the scale-dependence in Eq. ^ , we 
correct the Green function for the modified Newtonian 
model: 



yp,q,r ^ 



1 + a 
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Note that k in Eq. (PD|) needs to be given in the form, 
consistently with the Green function itself, as 



M 



k{p, q,r) ^ — { [sin^(7rp/A/) + sin^ (vrg/A/) 
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(41) 



We use the above Green function, evolve the system from 
z = 50 to 0, and make mock galaxy samples in the man- 
ner described in the next subsection. 



B. Observational data and mock samples 

For definiteness, we choose a volume-limited sample of 
SDSS galaxies whose r-band magnitude is in the range 
of (-21.0, -20.0) from those described in Hikage et al. 
[Hi. The redshift range is 0.044 < z < 0.103, the sur- 
vey volume, Kiamp, is 9.20 x 10^(/i~-^Mpc)^, and the total 
number of galaxies is 44,636. We made sure that using 
the other volume-limited samples with different magni- 
tude ranges [s^ does not significantly affect the results 
of our analysis below. 

We generate 24 mock catalogues from our A^-body sim- 
ulation data. The mock catalogues take into account var- 
ious observational effects such as survey geometry, the 
number density, and redshift distortion (peculiar veloc- 
ities of simulations particles are assigned to the mock 
galaxies ) [S^. In order to account for the effect of sur- 
vey geometry, we distribute random particles within the 
survey volume and correct for the boundary effect follow- 
ing the prescription of Feldman, Kaiser and Peacock [sj ■ 
We subtract fluctuations of the random particles which 
are within the survey volume, <5fc.random: 



^fc, random ■ 



(42) 



While this prescription is fairly empirical and may not 
completely account for the effect of the survey geom- 
etry, it yields a robust estimate at scales of our main 
interest here, k ~ 0.1/iMpc~^. When we calculate the 
power spectrum and bispectrum for SDSS galaxies and 
the mock catalogues, we use the above "corrected" den- 
sity, Sk. 



VI. CONSTRAINTS FROM POWER 
SPECTRUM 

We first compare the power spectra used the Peacock- 
Dodds prescription and those from numerical simula- 
tions. In Fig. [1] we plot the mass power spectra in 
real space (left panels) and in redshift space (right pan- 
els) . The predictions from perturbation theory agree well 
with the results of A^-body simulations.. Note that in the 
Newtonian case, the predicted power spectra with bi = 1 
are already in reasonable agreement with the observed 
power spectrum of SDSS galaxies. Our simulation re- 
sults are also consistent with those of Stabenau and Jain 

The panels on the right side in Fig. [1] show the power 
spectra of our mock "galaxies" . In each panel, the dot- 
ted line indicates the non-linear power spectrum in real 
space, which is the same in the corresponding left panel 
and shown for comparison. The redshift-space power 
spectrum of the SDSS volume- limited sample is shown by 
cross symbols. To include effects of redshift space distor- 
tion in our model, we use the formula derived of Magira, 
Jing and Suto [331 [equation (12) in their paper]. On lin- 
ear scales {k < O.l/iMpc""'^), the Kaiser effect is clearly 
seen as an enhanced power with respect to the real space 
power spectrum. It is worth mentioning that the plotted 
power spectra show substantial variations on the largest 
scales (k < 0.03ft-Mpc^^), which are presumably due to 
the somewhat complex survey geometry. 

To derive constraints on a, A using the calculated 
power spectra, we apply the Ax^ statistic. We treat 
the linear bias parameter bi as a free parameter in or- 
der to adjust the overall amplitudes of the power spectra 
between the predictions/simulations and the SDSS data. 
This normalization allows us to use the shape of the power 
spectra to detect possible deviations from the Newtonian 
case. 

We calculate as 



X 



E 



[Pjh) ~ Psussih)]' 
<y^{ki) 



(43) 



where PsDSs(fci) is the SDSS galaxy power spectrum. We 
use the predicted power spectra P{ki) and the variance 
of the SDSS data, a'^{ki), to calculate i^i real space, 
while for the same analysis in redshift space, we use those 
power spectra with the variance of mock galaxy samples 
to represent the cosmic variance in redshift space. 

We compute the relative confidence level of a and A 
with respect to their best-fit values assuming that 

Ax^ia, A) EE X^(a, A,6*,local mm) - X^(amin, Amin,6*,min) 

(44) 

follows the ^^distribution for 2 degrees of freedom. In Eq. 
(|44)) . aniin, Amin and 5*_niin dcuotc their best-fit values 
which globally minimize the value of x^, while 6*jocai min 
is the value that minimizes the x^ for a given set of values 
of a and A. 
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FIG. 1: The panels of tlie left side show the power spectra in real space. The adopted model parameters are (a) Newtonian 
(a = 0.0), (b)a — +1.0 and A — 5/i~^Mpc, (c) a = —1.0 and A = 5/i~^Mpc. Dots with vertical and horizontal error-bars are 
the power spectrum of SDSS galaxies from Tegmark et al. 01 • Dots with only vertical error-bars indicate results of N-body 
simulations. Dotted and solid lines are linear and non-linear power spectrum, respectively. We denote by k, the length scale 
of the mean inter-particle separation in our simulations, which is given by k, = 0.5 ■ 2-k ■ A^p^'^/Lbox- The simulation results 
are reliable at fc < In the right panels, we plot the power spectra for mock "galaxies" generated from our simulation. The 
parameters for (d), (e), (f) are the same as for (a), (b), (c), respectively. The dotted and solid lines in the right panels are 
non-linear power spectrum in real space (, which are the same as solid lines in the left panels) and redshift space, means the 
one- dimensional velocity dispersion calculated from simulation data. Cross symbols are the power spectra of the volume-limited 
sample of SDSS galaxies 



Figure shows the contours of Ax^ia, A). The results 
from TV-body simulations in real space are shown in panel 
(a) and (b). These differ only in the simulation box size, 
500 /i-^Mpc for (a) and 1000 h-^Mpc for (b). Hence the 



range of k used to derive constraints is slightly different, 
We also show the result from the real-space Peacock- 
Dodds prediction by thin dotted lines using the the same 
range of k consistently with the simulations. Clearly, the 
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FIG. 2: Constraints on a and A obtained using (a) simula- 
tions with I/box = 500/!.~^Mpc, (b) simulations with Lbox ~ 
1000/i^^Mpc, and (c) mock galaxy catalogues. The range of 
k used to compute is indicated in each panel. Thick solid, 
dotted, thick dotted lines indicate the limits at la, 2a and 3a 
confidence levels. For comparison, we plot the corresponding 
la, 2a and 3a confidence levels using the real-space Peacock- 
Dodds prediction in thin dotted lines. 



results of the perturbation theory and that of our numer- 
ical simulations are consistent with each other, putting 
quite similar constraints on a and A. 

The bottom panel (c) in Figure shows the constraints 
from our mock galaxy samples in redshift space. The 
range of k used in the analysis is 0.03 < k < 0.20/iMpc~^. 



The constraint is slightly less tight than those from 
perturbation theory and iV-body simulations. This is 
mainly because we discard the data points at large scales 
k ^ 0.01/iMpc~^ where the deviations from the Newto- 
nian case are most significant. Nevertheless models with 
I a I > 1 are still excluded at a 2-3(7 confidence level for 
A ~ 10ft-~^Mpc. For reference, we also plot the con- 
tours based on the real-space Peacock-Dodds prediction 
by thin dotted lines. 



VII. CONSTRAINTS FROM BISPECTRUM 

We further derive constraints on the modified New- 
tonian model extending the analysis to the three-point 
statistics. Specifically we use (conventional) bispectrum, 
B{ki,k2), defined in Eq. ([5^ . and reduced bispectra Q 
and p*^^^ defined as 



Q(fci,fc2) 



B{ki,k2) 



P(fcl)P(fc2)+P(fc2)P(fc3) 



and 



B{ki,k2) 



Vl4ampP(fcl)P(fc2)-P(fc3) 



P(fc3)F(fci)' 
(45) 



(46) 



where k^ = — fei — fe2, ki = \ki\, and K;amp is the sam- 
pling volume. The latter quantity p^^-* is the probabil- 
ity density function of phase sum for a density field, 
Oki+Ok^ +^fc3 [5k = |5fe|ex2(i6'fc)], studied in Matsubara 
poj and Hikage et al. [30, |41|. In this paper, we consider 
only isosceles triangles in fc-space that satisfy the relation 
k = ki = k2 with angle ip defined as 



ki ■ k2 

kik2 



(47) 



In the following analysis, we usep^'^^ to give constraints 
on the deviation from Newtonian gravity. This is because 
p('^) consists only of Fourier-phase informations and thus 
their constraints have good complementarity with those 
from P{k), which is defined as the square of the Fourier 
amplitudes. 



A. Linear bias model with b2 — 

Let us consider first linear bias model [62 = in 
Eq. (21)]. Figure [3] plots the bispectra B,Q,p^^^ in real 
space (left panels) and in redshift space (right panels) 
for Lbox =500 h~^Mpc simulations. The survey volume 
is set to be (ibox)'^ in Eq- (011) • The bispectra at small 
if are dominated by various nonlinear effects, whereas 
there are substantial uncertainties at large ip because 
of the small number of Fourier modes sampled. Given 
those, the agreement between predictions from perturba- 
tion theory (dashed lines) and iV-body simulation data 
(solid circles with error-bar) is very satisfactory. 
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FIG. 3: Bispectra B{k), Q{k), and p'''^(fc) from top to bottom as a function of ip measured in real space (left panels) and in 
redshift space (right panels). The range of fc = |fei| = |fe2| is indicated in each panel. The value of A is fixed as 5/i~^Mpc. The 
dashed, solid, dotted lines show the perturbation predictions in real space for a — +1.0, 0.0 (Newtonian) and —1.0, respectively. 
Symbols with error-bars show the results of simulations. Open circle, filled circle, cross symbol also mean a = +1.0, O.O(Newton) 
and —1.0, respectively. 



The right panels of Fig. [3] shows the bispectra in red- 
shift space. There, the results from our mock samples 
are shown by symbols with error bars. For comparison, 
we also show the results from perturbation theory in real 
space. In Fig. [2Id), Kaiser effect is clearly seen as a 
enhance at small ip. 

We further examine the dependence of the bispectra on 
A. Figure [5] compares the bispectra for different values of 



A. We have set a — 0.5 (left panels) and a ~ —0.5 (right 
panels). 

Figure [3] shows p^^^^ (fc) for the volume-limited SDSS 
catalogue and for our mock samples at fc in the range of 
0.072/iMpc"^ < fc < 0.103/iMpc"\ They have a very 
similar shape, but their amplitude depends systemati- 
cally on the value of a, the degree of deviations from the 
Newtonian case. 
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FIG. 4: Bispectra B{k), Q{k), and from top to bottom as a function of measured in real space; right: a — 0.5, left: 

a = —0.5. The dotted, dashed, and solid lines show the perturbation predictions in real space for A = 5/i^^Mpc, A = 10/i^^Mpc, 
and A = cxD(Newtonian), respectively, while symbols indicate the corresponding simulation results. 



Figure [S] plots constraints on the (a, A) plane derived 
from the Ax^ fit to the SDSS bispectrum using p'^'^^ and 
assuming a linear bias (62 = 0). The constraints from the 
bispectrum arc fairly consistent with, but sHghtly more 
stringent than, those from the power spectrum, which in- 
dicates the complementary role of the higher-order clus- 
tering statistics. 



B. The effect of non-linear biasing 



In reality, however, it may be more appropriate to ana- 
lyze the higher-order clustering statistics adopting a non- 
linear bias model. In the case of the bispectrum, it im- 
plies to introduce the quadratic biasing parameter 62 [see 

Eq. (121)]. In this bias model, the relation of pg for galax- 
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FIG. 5: Bispectra p'^' for SDSS galaxies (crosses), mock 
galaxy samples (solid circles with error bars), and perturba- 
tion theory predictions in real space (dotted line): (a) a = 
+ 1.0 and A = 5/i~^Mpc, (b) a = 0.0 (Newtonian), (c) 
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ies and p^^'' for mass reduces to 



(3)/ 



bi 



fiPl,P2,P3 



(48) 



P1P2 + P2P3 + P3P1 , ^ 

f{Pl,P2,P3) ^ ' ' p J ^, (49) 

where Pi = P{h) for i = 1, 2, 3 [s^. 

Previous papers [s^, suggest that a simple linear 



FIG. 6: Constraints on the a — A plane from the p^^^ analysis 
assuming 62 = 0. The range of k is from 0.07 to 0.20 /iMpc~^. 
Solid, dashed, thick dotted lines indicate la, 2a and 3a" con- 
fidence levels. Thin dotted lines are the same as those in Fig. 
ETa). 
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bias model in the Newtonian gravity model describes 
well the clustering of the volume-limited sample of SDSS 
galaxies, i.e., 62 ~ and 61 ~ 1. We now repeat the 
similar analysis in the modified Newtonian model. 

Figure[7Ja) indicates constraints on the (a. A) plane by 
treating 62/&1 as a free parameter, which should be com- 
pared with Figure [S] for 62 = 0. The regions below the 
contours are excluded with the corresponding confidence 
level. Naturally the bispectrum alone does not constrain 
(a. A) significantly in this generalized model. While the 
a — models are excluded with a la confidence level, 
the conclusion is not statistically significant. In turn, 
however, we can derive constraints on the value of 62/^1 
for the modified gravity model by combining the con- 
straints from power spectrum (independent of the value 
of 1)2 1 hi). Figure [Tfb) shows the contours of the best-fit 
value of 62/fci that gives the minimum for p^'^\k) on 
the plane. Figure [Tjb) suggests that 62/&1 should satisfy 
—0.4 < 62/61 < 0.3, which is the first constraint on the 
quadratic biasing parameter in the modified Newtonian 
model. 



This method allows us to compute the clustering statis- 
tics in redshift space and taking account of various ob- 
servational effects such as survey geometry as well. The 
resulting constraints from power spectrum are consistent 
with those obtained in our earlier work, indicating the va- 
lidity of the previous empirical modeling of gravitational 
nonlinearity in the modified Newtonian model. If linear 
biasing is adopted, the bispectrum of the SDSS galaxies 
yields constraints very similar to those from the power 
spectrum. If we allow for the nonlinear biasing instead, 
we find that the ratio of the quadratic to linear biasing 
coefficients, 62/61, should satisfy —0.4 < 62/61 < 0.3. in 
the modified Newtonian model. 

Future observations will exploit large ground-based 
telescopes to probe the matter density distribution by 
weak gravitational lensing. Combined with data from 
galaxy redshift surveys, lensing observations will provide 
invaluable informations on galaxy bias. Then it will be 
possible to put more stringent constraints on deviations 
from Newton's law of gravity at cosmological scales, us- 
ing the methodology presented in the this paper. 



VIII. SUMMARY 

We have derived constraints on possible deviations 
from Newtonian gravity using the power spectrum and 
the bispectrum of Sloan Digital Sky Survey galaxies. Our 
model assumes an additional Yukawa-like term with two 
parameters that characterize the amplitude, a, and the 
length scale. A, of the modified gravity. We have pre- 
dicted the power spectrum and the bispectrum using 
two different methods, the perturbation theory and di- 
rect A^-body simulations, and found the good agreement 
in real space as long as the biasing between galaxies and 
mass is neglected. In order to take the biasing effect into 
consideration, we adopt a quadratic biasing model. By 
comparing with the mock catalogues constructed from 
our simulations, we have derived constraints on a and A. 
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